Machine Learning Lab 04: Multinomial Logistic Regression¶

Generalized Linear Models¶

From lecture, we know that members of the exponential family distributions can be written in the form $$p(y;\eta) = b(y)e^{(\eta^\top T(y)-a(\eta))},$$ where

• $\eta$ is the natural parameter or canonical paramter of the distribution,
• $T(y)$ is the sufficient statistic (we normally use $T(y) = y$),
• $b(y)$ is an arbitrary scalar function of y, and
• $a(\eta)$ is the log partition function. We use $e^{a(\eta)}$ just to normalize the distribution to have a sum or integral of 1.

Each choice of $T$, $a$, and $b$ defines a family (set) of distributions parameterized by $\eta$.

If we can write $p(y \mid \mathbf{x} ; \theta)$ as a member of the exponential family of distributions with parameters $\mathbf{\eta}$ with $\eta_i = \theta^\top_i \mathbf{x}$, we obtain a generalized linear model that can be optimized using the maximum likelihood principle.

The GLM for the Gaussian distribution with natural parameter $\eta$ being the mean of the Gaussian gives us ordinary linear regression.

The Bernoulli distribution with parameter $\phi$ can be written as an exponential distribution with natural parmeter $\eta = \log \frac{\phi}{1-\phi}$. The GLM for this distribution is logistic regression.

When we write the multinomial distribution with paremeters $\phi_i > 0$ for classes $i \in 1..K$ with the constraint that $$\sum_{i=1}^{K} \phi_i = 1$$ as a member of the exponential family, the resulting GLM is called multinomial logistic regression. The parameters $\phi_1, \ldots, \phi_K$ are written in terms of $\theta$ as $$\phi_i = p(y = i \mid \mathbf{x}; \theta) = \frac{e^{\theta^\top_i \mathbf{x}}}{\sum_{j=1}^{K}e^{\theta^\top_j \mathbf{x}}}. $$

Optimizing a Multinomial Regression Model¶

In multinomial regression, we have

1. Data are pairs $\mathbf{x}^{(i)}, y^{(i)}$ with $\mathbf{x}^{(i)} \in \mathbb{R}^n$ and $y \in 1..K$.

2. The hypothesis is a vector-valued function $$\mathbf{h}_\theta(\mathbf{x}) = \begin{bmatrix} p(y = 1 \mid \mathbf{x} ; \theta) \

                                     p(y = 2 \mid \mathbf{x} ; \theta) \\
                                     \vdots \\
                                     p(y = K \mid \mathbf{x} ; \theta) \end{bmatrix},$$
   

   where $$p(y = i \mid \mathbf{x}) = \phi_i = p(y = i \mid \mathbf{x}; \theta) = \frac{e^{\theta^\top_i \mathbf{x}}}{\sum_{j=1}^{K}e^{\theta^\top_j \mathbf{x}}}. $$

We need a cost function and a way to minimize that cost function. As usual, we try to find the parameters maximizing the likelihood or log likelihood function, or equivalently, minimizing the negative log likelihood function:

where $$\begin{eqnarray} J(\theta) & = & - \ell(\theta) \\ & = & - \sum_{i=1}^m \log p(y^{(i)} \mid \textbf{x}^{(i)} ; \theta). \end{eqnarray}$$

Now that we know what is $J(\theta)$, let's try to find its minimimum by taking the derivatives with respect to an arbitrary parameter $\theta_{kl}$, the $l$-th element of the parameter vector $\theta_k$ for class $k$. Before we start, let's define a variable $a_k$ as the linear activation for class $k$ in the softmax function: $$ a_k = \theta_k^\top \mathbf{x}^{(i)}, $$ and rewrite the softmax more conveniently as $$ \phi_k = \frac{e^{a_k}}{\sum_{j=1}^K e^{a_j}}. $$ That makes it a little easier to compute the gradient: $$\begin{eqnarray} \frac{\partial J}{\partial \theta_{kl}} & = & - \sum_{i=1}^m \frac{1}{\phi_{y^{(i)}}} \frac{\partial \phi_{y^{(i)}}}{\partial \theta_{kl}}. \\ \end{eqnarray}$$ Using the chain rule, we have $$\frac{\partial \phi_{y^{(i)}}}{\partial \theta_{kl}} = \sum_{j=1}^K \frac{\partial \phi_{y^{(i)}}}{\partial a_j} \frac{\partial a_j}{\partial \theta_{kl}}$$ The second factor is easy: $$ \frac{\partial a_j}{\partial \theta_{kl}} = \delta(k=j)x^{(i)}_l. $$ For the first factor, we have $$\begin{eqnarray} \frac{\partial \phi_{y^{(i)}}}{\partial a_j} & = & \frac{ \left[ \delta(y^{(i)}=j)e^{a_j} \sum_{c=1}^K e^{a_c} \right] - e^{a_j} e^{a_j} }{\left[ \sum_{c=1}^K e^{a_c} \right]^2} \\ & = & \delta(y^{(i)}=j) \phi_j - \phi_j^2 \end{eqnarray}$$

Substituting what we've derived into the definition above, we obtain $$ \frac{\partial J}{\theta_{kl}} = - \sum_{i=1}^m \sum_{j=1}^K (\delta(y^{(i)}=j) - \phi_j) \frac{\partial a_j}{\partial \theta_{kl}}. $$

There are two ways to do the calculation. In deep neural networks with multinomial outputs, we want to first calculate the $\frac{\partial J}{\partial a_j}$ terms then use them to calculate $\frac{\partial J}{\partial \theta_{kl}}$.

However, if we only have the "single layer" model described up till now, we note that $$\frac{\partial a_j}{\partial \theta_{kl}} = \delta(j=k) x^{(i)}_l,$$ so we can simplify as follows: $$\begin{eqnarray} \frac{\partial J}{\theta_{kl}} & = & - \sum_{i=1}^m \sum_{j=1}^K (\delta(y^{(i)}=j) - \phi_j) \frac{\partial a_j}{\partial \theta_{kl}} \\ & = & - \sum_{i=1}^m \sum_{j=1}^K (\delta(y^{(i)}=j) - \phi_j) \delta(j=k) x^{(i)}_l \\ & = & - \sum_{i=1}^m (\delta(y^{(i)}=k) - \phi_k) x^{(i)}_l \\ \end{eqnarray}$$

Put It Together¶

OK! Now we have all 4 criteria for our multinomial regression model:

1. Data are pairs $\mathbf{x}^{(i)}, y^{(i)}$ with $\mathbf{x}^{(i)} \in \mathbb{R}^n$ and $y \in 1..K$.

2. The hypothesis is a vector-valued function $$\mathbf{h}_\theta(\mathbf{x}) = \begin{bmatrix} p(y = 1 \mid \mathbf{x} ; \theta) \

                                     p(y = 2 \mid \mathbf{x} ; \theta) \\
                                     \vdots \\
                                     p(y = K \mid \mathbf{x} ; \theta) \end{bmatrix},$$
   

   where $$p(y = i \mid \mathbf{x}) = \phi_i = p(y = i \mid \mathbf{x}; \theta) = \frac{e^{\theta^\top_i \mathbf{x}}}{\sum_{j=1}^{K}e^{\theta^\top_j \mathbf{x}}}. $$

3. The cost function is $$J(\theta) = - \sum_{i=1}^m \log p(y^{(i)} \mid \textbf{x}^{(i)})$$

4. The optimization algorithm is gradient descent on $J(\theta)$ with the update rule $$\theta_{kl}^{(n+1)} \leftarrow \theta_{kl}^{(n)} - \alpha \sum_{i=1}^m (\delta(y^{(i)}=k) - \phi_k) x^{(i)}_l.$$

Multinomial Regression Example¶

The following example of multinomial logistic regression is from Kaggle.

The data set is the famous Iris dataset from the UCI machine learning repository.

The data contain 50 samples from each of three classes. Each class refers to a particular species of the iris plant. The data include four independent variables:

1. Sepal length in cm
2. Sepal width in cm
3. Petal length in cm
4. Petal width in cm

The target takes on one of three classes:

1. Iris Setosa
2. Iris Versicolour
3. Iris Virginica

To predict the target value, we use multinomial logistic regression for $k=3$ classes i.e. $y \in \{ 1, 2, 3 \}$.

Given $\mathbf{x}$, we would like to predict a probability distribution over the three outcomes for $y$, i.e., $\phi_1 = p(y=1 \mid \mathbf{x})$, $\phi_2 = p(y=2 \mid \mathbf{x})$, and $\phi_3 = p(y=3 \mid \mathbf{x})$.

The phi function returns $\phi_i$ for input patterns $\mathtt{X}$ and parameters $\theta$.

Tips for using intellisense: Shift+TAB

The grad_cost function gives the gradient of the cost for data $\mathtt{X}, \mathbf{y}$ for class $j\in 1..k$.

Exercise 1.1 (5 points)¶

Create a function to load data from Iris.csv using the Pandas library and extract y from the data.

You can use the Pandas 10 minute guide to learn how to use pandas.

Expected result: \ SepalLengthCm SepalWidthCm PetalLengthCm PetalWidthCm \ 0 5.1 3.5 1.4 0.2\ 1 4.9 3.0 1.4 0.2\ 2 4.7 3.2 1.3 0.2\ 3 4.6 3.1 1.5 0.2\ 4 5.0 3.6 1.4 0.2\ 0 Iris-setosa\ 1 Iris-setosa\ 2 Iris-setosa\ 3 Iris-setosa\ 4 Iris-setosa\ Name: Species, dtype: object

Exercise 1.2 (10 points)¶

Partition data into training and test sets

• No need to use random.seed function!
• Ensure that the train set is 70% and the test set is 30% of the data.
• Encode the labels in the y attribute to be integers in the range 0..k-1.

Expected result: (*or similar*)\ X_train.shape (105, 4)\ X_test.shape (45, 4)\ y_train.shape (105,)\ y_test.shape (45,)\ y_labels_name: ['Iris-setosa' 'Iris-versicolor' 'Iris-virginica'] \ y_labels_new: [0, 1, 2]

SepalLengthCm SepalWidthCm PetalLengthCm PetalWidthCm\ 0 6.4 2.8 5.6 2.2\ 1 6.7 3.3 5.7 2.1\ 2 4.6 3.4 1.4 0.3\ 3 5.1 3.8 1.5 0.3\ 4 5.0 2.3 3.3 1.0\ Species\ 0 2\ 1 2\ 2 0\ 3 0\ 4 1

Exercise 1.3 (5 points)¶

Train your classification model using the gradient_descent function already provided. You might also play around with the gradient descent function to see if you can speed it up!

Expected result: (*or similar*)\ [[ 1.17632192 1.32360047 1.83204165 -0.20224445 0.44039155]\ [ 1.10140069 1.13537321 0.74833178 1.21907866 0.82567377]\ [ 0.72227738 0.54102632 0.41962657 1.98316579 1.73393467]]\ \ (3, 5)

Exercise 1.4 (5 points)¶

Let's get your model to make predictions on the test data.

Expected result: (*or similar*)\ 45 \ [2, 0, 2, 0, 0, 0, 2] \

<class 'int'>

Exercise 1.5 (5 points)¶

Estimate accuracy of model on test data

Expected result: should be at least 0.8!

On your own in lab¶

We will do the following in lab:

1. Write a function to obtain the cost for particular $\mathtt{X}$, $\mathbf{y}$, and $\theta$.
2. Plot the training set and test cost as training goes on and find the best value for the number of iterations and learning rate.
3. Make 2D scatter plots showing the predicted and actual class of each item in the training set, plotting two features at a time. Comment on the cause of the errors you observe. If you obtain perfect test set accuracy, re-run the train/test split and rerun the optimization until you observe some mistaken predictions on the test set.

Exercise 2.1 (15 points)¶

1. Write a function to obtain the cost for particular $\mathtt{X}$, $\mathbf{y}$, and $\theta$. Name your function my_J() and implement

1. Implement my_grad_cost using your my_J function

Expect result: (*or similar*)\ intercept 0.009524\ SepalLengthCm 0.316825\ SepalWidthCm -0.091429\ PetalLengthCm 0.780000\ PetalWidthCm 0.329524\ dtype: float64\ 37.352817814715735

1. Implement my_gradient_descent using your my_grad_cost function

Expected result: (*or similar*)\ [[1.00001186 0.99618853 1.00183642 0.9889817 0.99528923]\ [1.00009697 1.0011823 0.99883395 1.00316763 1.00083055]\ [0.99987915 1.00255606 0.99929351 1.00779768 1.00386218]]\ [114.00099216453735, 113.89036233839263, 113.78163144339288, 113.67472269747496, 113.56956268162737]\ 37.352817814715735

Exercise 2.2 (20 points)¶

1. Plot the training set and test cost as training goes on and find the best value for the number of iterations and learning rate.
2. Make 2D scatter plots showing the predicted and actual class of each item in the training set, plotting two features at a time. Comment on the cause of the errors you observe. If you obtain perfect test set accuracy, re-run the train/test split and rerun the optimization until you observe some mistaken predictions on the test set.

From the plots, it is seen that the learning rate of 0.1 and 450 number of iternations produced the best results. The cost reduced relatively more quickly in this configuration than other combinations. However, at one point during the training, the cost increased and performed worse than the model trained with a lower learning rate of 0.07 as observed from the graph. Though it started decreasing again shortly, the cost did not decrease smoothly. It is only after a certain iternation, did the cost of model decreased lower than the cost of other model and also smoothly. From this, it can be concluded that increasing learning more than this point might cause the model learn from data non-smoothly and might perform worse for a certain amount of iteration.

It is noted that the cost for model which performed the best would still decrese a bit and is still yet to reach its convergence. Maybe training for a longer number of iternations would slightly improve the model. One optimal algorithm which can be used for this problem is early stopping. It is basically a method that stops the training once model stops improving to avoid overfitting. Using this method, larger number of iternations can be defined and hence eliminates the task for defining optimal number of iteratinons to large extend.

Expected result: (*Yours doesn't have to be the same!*)

Take Home Solution¶

Reference¶

https://www.kaggle.com/parulpandey/palmer-archipelago-antarctica-penguin-data

The dataset which I used for take-home exercise was taken from 'Palmer Archipelago (Antarctica) penguin data' in Kaggle. It is a great dataset for practising data exploration & visualization and is seen as an iternative to iris dataset. For my exeriment, based on the information provided on penguins, I tried to create a suitable model which will guess the species of the penguin. The species given in this dataset are chinstrap, gentoo, and adelie.

Note: There were two datasets provided from 'Palmer Archipelago (Antarctica) penguin data'. One was the simplified version of the other data. The one I used was the simplified version since it contained all the important information and omits unnecessary details.

Loading Data¶

Handling Missing Values¶

Handling Missing Value for Sex Category¶

Handling Missing Value for culmen_length_mm Category¶

Handling Missing Value for Other Categories in a Similar Way¶

Converting Categorical Columns into Numeric Entries¶

Exploring 'sex' column¶

Exploring 'island' column¶

Normalizing Data¶

Data for columns 'culmen_length_mm','culmen_depth_mm','flipper_length_mm', and 'body_mass_g' were normalized since taking exponential of the dot product of training data and theta results in 'inf' due to their large values.

Splitting Data into Train and Test Dataset¶

Training the Model¶

Plot¶

Summary¶

1) Data from 'Palmer Archipelago (Antarctica) penguin data' was uploaded

2) Missing Values in the data were filled in

3) Data for columns 'culmen_length_mm','culmen_depth_mm','flipper_length_mm', and 'body_mass_g' were normalized

4) Categorical columns were converted into numeric entries

5) Data was split into train and test datasets

6) Model was trained with 4 different learning rates

Conclusion and Remarks¶

From iternation vs. cost plot, it is observed that increasing learning rate allows the model to learn from the data more quickly for a few number of iterations and hence reaches its convergence sooner. However, as the training progessed, the model trained with lower learning rate seemed to be catching up with the model with higher learning rate. From this observation, one way which can be used to train the model more effeciently is to first begin with a relatively large learning rate. Then as the training progresses, the learning rate of the model decreases gradually to increase the probablity of the model converging in local minima of cost function. Early stopping, discusssed briefly in exercsie 2.2, can also be added to further improve the training so longer number of iternations can be used.